rm(list=ls(all=TRUE))
library(stats)
# Settings: n=10,000, fraud=.5, districts(j)=100 
# Calibration: alphaA=400, alphaB=320, stolen(gamma)=.3, used(delta)=1.2
simulations=function(n,fraud,districts,alphaA,alphaB,stolen,used){

data=matrix(NA,n, 9,dimnames=list(seq(1,n), c("fraud","votesA","votesB", "newA", "newB","meanAf","meanBf","Af1","Bf1")))

#Proportion of "fraudulent" simulations
dummie=rbinom(n,1,fraud)
data[,1]=dummie

for (j in 1:n)
{
#Generating simulated vote countes (taking out negative values)
# Benford Variate
XA=10^runif(districts,0,1)
XB=10^runif(districts,0,1)
# Votes for party i (equation 2)
votesA=as.integer(alphaA*XA)
votesB=as.integer(alphaB*XB)

#stolen votes
lostvotesA=as.integer(stolen*votesA)
#filter
lostAfilter=lostvotesA*dummie[j]

#New votes for A
newvotesA=votesA-lostAfilter

#New counts after the "fraud"
##1. All the votes lost by party A go to party B
newvotesB=as.integer(votesB+(used*lostAfilter))

# Generate first digits of data
firstA=as.numeric(substring(newvotesA,1,1))
firstB=as.numeric(substring(newvotesB,1,1))

data[j,2]=(sum(votesA))/districts
data[j,3]=(sum(votesB))/districts
data[j,4]=(sum(newvotesA))/districts
data[j,5]=(sum(newvotesB))/districts
data[j,6]=mean(firstA,na.rm=TRUE)
data[j,7]=mean(firstB,na.rm=TRUE)
data[j,8]=(sum(firstA==1,na.rm=TRUE))/districts
data[j,9]=(sum(firstB==1 ,na.rm=TRUE))/districts
}

print(data)
write.csv(data,file="/Users/sebastian/Mac/Argentina/Fraud/vote_simulation",col.names=TRUE)
}
